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Stable fluid and solid particle phases are essential to the simulation of continuum fluids and solids 
using Smooth Particle Applied Mechanics. We show that density-dependent potentials, such as 
= 2 XX/ 9 — Po) 2 ) along with their corresponding constitutive relations, provide a simple means 
for characterizing fluids and that a special stabilization potential, <l?vp = ^XX^/o) 2 , not only 
stabilizes crystalline solid phases (or meshes) but also provides a surface tension which is missing 
in the usual density-dependent-potential approach. We illustrate these ideas for two-dimensional 
square, triangular, and hexagonal lattices. 
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I. SMOOTH PARTICLE APPLIED MECHANICS 

Smooth particle applied mechanics — SPAM — was discovered about thirty years ago 
It has become a useful tool in simulating gases, fluids, and solids, and holds particular 
promise for problems involving large high-speed deformation and failure. The main advan- 
tage of the method is simplicity. SPAM closely resembles atomistic molecular dynamics and, 
in a variety of cases[3j], the SPAM particle trajectories are isomorphic to those of molecular 
dynamics. 

A smooth-particle code is less-complicated than typical grid-based continuum codes be- 
cause the smooth-particle method evaluates spatial gradients in a particularly simple way, 
explained in more detail below. The main disadvantages of the method are instability in 
tension j^J and the lack of surface tension^. The present work introduces an idea — density- 
gradient potentials — designed to address those problems. 

The basic smooth-particle approach is to represent all continuum properties (the density 
p, the velocity v, the stress tensor cr, ...) as interpolated sums of particle properties, where 
the particles are described by "weight functions", expressing the range of influence of the 
particles in space. The simplest weight function satisfying five desirable conditions — (i) 
normalization, (ii) finite range h, (iii) a maximum at the origin, and (iv and v) two continuous 
derivatives everywhere — is Lucy's. Normalized for applications in two-dimensional space 
Lucy's weight function is as follows: 
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27rrw(r)dr = 1 . 

The density at any point r is defined as the sum of all the particle contributions at that 
point: 

p(r) = '^2m,jw(r — rj) , 

j 

so that the density associated with Particle i is 

j 

Other continuum properties at location r are likewise calculated as sums over nearby parti- 
cles: 

f(r)p(r) = Y, m jfjw(\r - r,-|) — ► 



/( r ) = [f(r)p(r)]/ p(r) = 
j j 

It is important to note that the interpolated function f(r) at the location of Particle j is 
typically different to the particle property fj at that point. 

Beyond using the point properties {fj} to define the field properties f(r) this approach 
has the crucial advantage that gradients of the field properties translate into simple sums of 
particle quantities: 

V(p/)r = (/V P ) r + ( P V/) r = 



Y / m j f j w(\r - rj \) 



J2 m jfj^rw(\r - rj\) . 



Expressions for the gradients of density, velocity, stress, and energy make it possible to 
express the partial differential equations of continuum mechanics as ordinary differential 
equations for the evolution of the particle coordinates, velocities, stresses, and energies| 
The resulting "equation of motion" for the particles is 



[{mP/p 2 )i + {mP/p 



where Pj is the pressure tensor associated with Particle %. Where the pressure is hydrostatic, 
and slowly varying in space, it is noteworthy that the smooth-particle equations of motion 
are exactly the same as the equations of molecular dynamics, with the weight function w(r) 
playing the role of a pair potential. In the simple case that the internal energy depends only 
on volume (and not on temperature) the pressure is simply related to the internal energy 
per unit mass e: 

P = p 2 de/dp . 

Specifying the density dependence of either the pressure or the internal energy, in such a 
case, corresponds to giving a full description of the equilibrium equation of state. 



II. CONVERGENCE OF SMOOTH-PARTICLE AVERAGES 

The evolving time-and-space dependent smooth-particle sums converge to continuum 
mechanics as a many-particle limit, just as do the more usual grid-based approximations. 
The range of the smooth-particle weight function, h, corresponds to a few grid spacings. In 
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FIG. 1: Summed-up densities evaluated at regular-lattice particle sites for (from top to bottom at 
the left side of the plot) the square, triangular, and hexagonal lattices. The range of Lucy's weight 
function varies from 2 to 5 where the overall density of the lattice is unity and all the particles 
have unit mass. 

practice, for an error level of order one percent, the weight-function sums must include a few 
dozen particles. Fig. 1 shows the dependence of density sums, pi = J2j mjWij on the range 
of the weight function for three regular two-dimensional lattices. In typical applications, 



The density curves in the Figure correspond to an actual density of unity. The many 
crossings of the curves suggest that energetic flows of highly inhomogeneous fluids would 
exhibit a complex structure without any definite lattice structure while very slow flows 
might "freeze" into a least-energy crystalline form. Simulations of the Rayleigh-Benard 
problemfconvection driven by a temperature gradient in the presence of gravity) support this 
surmise p]. Low-energy, high-pressure simulations can actually "freeze", with the smooth 
particles forming a locked lattice structure rather than flowing. For fluids this freezing 
behavior is undesirable. In the next Section we consider the stress-free mechanical stability 
of the three simplest two-dimensional lattice structures. 



with h ~ 3JV/N the density errors are of order one percent. 
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III. PHASE INSTABILITY FROM DENSITY POTENTIALS 



When discrete particles are involved there can be difficulties in representing the smooth 
and continuous nature of fluid flows. By analogy with atomistic molecular dynamics, one 
would expect that regular lattice arrangements of particles would resist shear. In the atom- 
istic case in two space dimensions the shear modulus G is of the same order as the one-particle 
Hooke's-law force constant evaluated from the Einstein model: 

r ~ - d ^ 

" — ^Einstein — 



dx\ 

where x\ is the displacement of a single test particle, Particle 1, from its lattice site, with 
all the other particles fixed. For a sufficiently simple density-dependent potential we can 
estimate the one-particle force constant ^Einstein analytically. 
Let us illustrate for the simplest possible density potential, 

$ p = H \\pj - Po] 2 ; Pj = Y. 771 ^ > 

3 Z i 

where $ p is the total potential energy of the system, p is the target density minimizing 
that energy, and all the particle masses are set equal to unity, rrij — 1. The first derivative, 



dxi 
vanishes for 

x\ = -> pi = pj = po . 
The second derivative can be estimated by replacing the particle sum with an integral: 
<9 2 $ p ^ ,o & fxw'\\ , 90 



dx\ Y Jo \ r J 7nh 4 ' 

Fig. 2 shows that this analytic result closely resembles the detailed lattice sums for all three 
regular two-dimensional lattices. 

Nevertheless, our detailed investigation of this particular choice of fluid model, 

E = E \{P - Po) 2 ^P = p\p - Po) , 

revealed that this expectation of a shear strength varying as h~ A is unfounded. Instead, 
regular lattices, with a stress-free, density-based potential corresponding to an athermal 
fluid constitutive relation, show no shear resistance whatever! 




FIG. 2: Comparison of the exact summed- up Einstein force constant ^Einstein with the approximate 
integrated estimate as a function of the range h of Lucy's weight function. The top-to-bottom 
ordering of the curves is [integrated > triangular > square > hexagonal] . 



Numerical investigation shows that the square, triangular, and hexagonal lattices, ar- 
ranged at the target density po, are all unstable to small displacements. This can be shown 
by using lattice dynamics, elastic theory, or molecular dynamics. In every case the regular 
lattices are unstable to a variety of shear modes. 

The perfect-crystal elastic constants^, 0] for this potential can be calculated by two 
chain-rule differentiations of the potential <£> with respect to the elastic strains: 

<9 2 $ 



_ <9 2 $ _ d 2 $ 
C\\V = ; C12V = — — 



C 44 V 



du x 



du,, 



xx'-'^yy 



^-xx r% 1 '-yy 1 ^xy ., . 

ox ay ay ox 

Here u(r) = (u x , u y ) represents an infinitesimal displacement from the perfect-lattice con- 
figuration. The resulting elastic constants take the form of lattice sums: 

\ 2 

CuV = E|EW/r)] < 



de 2 

^^xy 



duv du, h 



C12V = £ YXA^'Mh YWW/r))* 
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FIG. 3: Variation of the bulk modulus B = C\\ + C\2 = 2Cn = 2Cyi with the range of the weight 
function h for three two-dimensional lattices. The top-to-bottom ordering of the curves at h = 2.1 
is [hexagonal > square > triangular]. 

C 44 V = £ (Efa/(ti//r)]y 

'"ij \j -^ij Vij ) "^ij "^i "Ej i Uij Hi Uj • 

For the square, triangular, and hexagonal lattices it is evident, by symmetry, that C\\ and 
C12 are equal and that C44 vanishes. The nonvanishing elastic constant Cn = C12 is exactly 
half the bulk modulus B. The range- dependence B{h) is shown in Fig. 3 for all three lattice 
structures. 




IV. PHASE STABILITY FROM DENSITY GRADIENTS 

The results of the preceding section show that the smooth-particle fluid model (correctly) 
is able to flow under an infinitesimal shear stress. For solids shear resistance is required. A 
simple potential supporting shear strength minimizes the gradient of the density: 

$VpOc^(Vp)5 . 
3 

This potential is minimized for regular lattices, in which there can (by symmetry) be no 
density gradient at the particle sites. For systems with free surfaces — the details are not 



FIG. 4: Hexagonal-lattice particle trajectories without the density-gradient potential. Initially 
the particle displacements were chosen randomly, with zero sum and with an initial rms value of 
y/ (5r 2 ) = 0.02. The range of Lucy's weight function is h = 3 with the density and the particle 
mass both chosen equal to unity. The elapsed time (40,000 Fourth Order Runge-Kutta timesteps 
dt = 0.05) is about 80 Einstein vibrational periods (^Einstein — 0.05). <£ p = I J2(p ~ Po) 2 - 

considered here, but are elaborated in a forthcoming book[l(J] — this potential also provides 
a surface tension, eliminating the tendency of smooth particles to form string-like phases. 
Figs. 4 and 5 illustrate the stability of the hexagonal lattice in the absence, and in the 
presence, respectively of the density-gradient potential. In the one example detailed here 
(which is typical of many we have investigated, with various sizes, initial conditions, and 
crystal structures) the individual particle trajectories with and without the density-gradient 
potential are shown. Evidently, by choosing the proportionality constants wisely, these 
potentials can be tuned to reproduce desired flow stresses for solids modelled with SPAM. 
This approach avoids many of the difficulties involved in integrating the smooth-particle 
equations for the stress rates, ({&} — > {ex}). 
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FIG. 5: Hexagonal-lattice particle trajectories mi/i the density-gradient potential <3>vp = 
^XXVp) 2 . Initial conditions and length of the simulation are identical to those of Figure 4. 
<D = <D p + <D Vp = \ U(p ~ Po) 2 + (Vp) 2 ]. 

V. CONCLUSIONS 

Density-dependent potentials can be used to simulate the behavior of either fluids or 
solids from the standpoint of smooth-particle simulation. By introducing density-gradient 
potentials strength and surface tension can be introduced, providing a useful model for 
solids. We believe that this idea will prove fruitful in a wide variety of high-strain-rate 
applications of smooth-particle methods. 
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